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ABSTRACT 

We have performed detailed population synthesis on a large number (2 x 10 7 ) of binary sys- 
tems in order to investigate the properties of massive double degenerate binaries. We have 
included new important results in our input physics in order to obtain more reliable estimates 
of the merging timescales and relative formation rates. These improvements include refined 
treatment of the binding energy in a common envelope, helium star evolution and reduced 
kicks imparted to newborn black holes. The discovery and observations of GRB afterglows 
and the identification of host galaxies have allowed comparisons of theoretical distributions 
of merger sites with the observed distribution of afterglow positions relative to host galax- 
ies. To help investigate the physical nature of short- and long-duration y-ray bursts (GRBs), 
we compute the distances of merging neutron stars (NS) and/or black holes (BH) from the 
centers of their host galaxies, as predicted by their formation scenario combined with motion 
in galactic potentials. Furthermore, we estimate the formation rate and merging rate of these 
massive double degenerate binaries. The latter is very important for the prospects of detecting 
gravitational waves with LIGO/VIRGO. We find that the expected detection rate for LIGO II 
is ~ 850 yr~' for galactic field sources and that this rate is completely dominated by merging 
BHBH binaries. Even LIGO I may detect such an event (~ 0.25 yr~'). Our preferred model 
estimate the Galactic field NSNS merger rate to be ~ 1 .5 x 10~ 6 yr _I . For BHBH systems this 
model predicts a merger rate of ~ 9.7 x 10~ 6 yr _1 . Our studies also reveal an accumulating 
numerous population of very wide orbit BHBH systems which never merge (t » THubbie)- 
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1 INTRODUCTION 

It has been recognized for a long time that the time duration of 
GRBs is bimodal: the majority (75%) of the bursts have a long 
duration with a mean of ~ 20 sec., the rest have a much shorter 
duration with a mean of only ~ 0.2 sec. (Mazets et al. 1981 and 
Hurley et al. 1992). Both of these two classes of GRBs have sim- 
ilar isotropic spatial distributions, but differ with respect to spec- 
tral hardness, fluence and temporal pulse properties (e.g. Kouve- 
liotou et al. 1993; Lee & Petrosian 1997 and Norris et al. 2000). It 
is natural to suggest that two distinct types of progenitors are re- 
sponsible for the existence of the two observed distinct populations 
of GRBs. 

Since it has been established from measurements of GRB red- 
shifts that their distance scales are cosmological (Metzger et al. 
1997), GRB scenarios require an energy release output of roughly 
(Q. y /4n) x 10 54 erg, where Q. y is the beaming angle. According to 
the internal shock model for the production of the observed y-rays 
(e.g. Piran 2000, and references therein), the duration of a burst 
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is very likely to be a direct measure of the time interval during 
which the powering engine is active. In both of the main scenar- 
ios the basic ingredient is a black hole surrounded by an accretion 
disk of nuclear matter threaded by a strong magnetic field. A large 
disk viscosity causes the matter to spiral inwards on a timescale of 
minutes. The poloidal part of the magnetic field in the inner disk, 
spinning at millisecond periods, will then accelerate a small amount 
of material (less than 10~ 5 M G ) in the form of a narrow and highly 
relativistic jet (this is the so-called Blandford-Znajek mechanism, 
e.g. Lee, Wijers & Brown 2000). Alternatively, the jet of material 
is produced by annihilating neutrinos along the disk axis. These 
neutrinos result from the release of binding energy of the rapidly 
accreted (> 0.01 M Q ) material. In either mechanism, the relativistic 
jet becomes a fireball and ultimately produces the GRB. 

The most promising models to account for the long and 
short duration bursts involve the collapse of a massive star (e.g. 
Woosley 1993 and MacFadyen & Woosley 1999) and the merg- 
ing of compact objects caused by the in-spiral due to gravitational 
wave emission in a tight binary (e.g. Goodman 1986; Eichler et al. 
1989; Paczynski 1990 and Meszaros & Rees 1992), respectively. 
These associations are supported by simulations of outburst dura- 
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tions which are estimated to be ~ several tens of seconds for the rel- 
ativistic outflow generated by the collapse of a massive star (Mac- 
Fadyen & Woosley 1999), and < 1 sec (Lee & Kluzniak 1998; 
Ruffert & Janka 1999) for the neutrino-driven wind of merging 
compact objects, respectively. It has been demonstrated (Meszaros, 
Rees & Wijers 1998) that both models are able to produce the re- 
quired GRB energies. 

An important method to test the above scenarios is by deter- 
mining the locations and environments in which the GRBs must 
occur according to their respective models. Collapsars (collapsing 
massive stars) are expected to be found close to their place of birth 
as a result of the short lifetime (a few Myr) of their progenitors. 
Therefore, these GRBs should occur in dense, dusty star-forming 
regions - i.e. inside their host galaxies. This hypothesis is in agree- 
ment with the observations of long-duration GRB afterglows: there 
seems to be a correlation between GRB locations and the UV light 
of their host galaxies (Sahu et al. 1997; Fruchter et al. 1998; Kulka- 
rni et al. 1998, 1999; Bloom, Kulkarni & Djorgovski 2002). Fur- 
thermore, an underlying supernova component (a "bump") seems 
to be present in the light curves of GRB 980326 (Bloom et al. 1999) 
and GRB 970228 (Reichart 1999; Galama et al. 2000), and recent 
detections of an iron line in the afterglow of several GRBs also 
provides evidence for the presence of a dense matter in the vicin- 
ity of the bursts (e.g. Yoshida et al. 1999; Piro et al. 1999, 2000; 
Vietri et al. 1999). 

All of the above mentioned afterglow observations only orig- 
inate from long-duration GRBs. So far, no similar types of obser- 
vations have been possible for the short bursts. It is the hope that 
observations of the short bursts will be possible in the near future. 
An analysis of the expected afterglow properties of these bursts 
is, for example, given by Perna & Belczynski (2002), while the 
prospect of future observations of these afterglows are discussed by 
Panaitescu, Kumar & Narayan (2001). The main topic to be inves- 
tigated is the distribution of merging compact binaries with respect 
to their host galaxies. The reason is that the two supernova explo- 
sions in the progenitor binaries result in systemic recoil velocities 
of typically a few hundred km s -1 , and therefore the merging com- 
pact binaries have the potential to travel far from their birth places 
- if their merger timescale is relatively large. However, the latter 
issue depends on the final orbital separation of the binaries, after 
the second star collapses to form a supernova. 

To estimate the merger timescales (and thus the offset of the 
associated GRBs from their birth places) it is necessary to per- 
form a binary population synthesis study of their progenitors. Such 
studies have been performed recently in a number of papers: e.g. 
Bloom, Sigurdsson & Pols (1999); Belczynski, Kalogera and Bu- 
lik (2002); Perna & Belczynski (2002); Belczynski, Bulik and 
Kalogera (2002); Sipior & Sigurdsson (2002). The applied forma- 
tion scenarios, input physics and subsequent results of these studies 
are not conclusive. Our major concern, however, is that all of the 
above papers have assumed rather simplified calculations e.g. for 
the very important common envelope and spiral-in phase by using 
a constant value for the so-called /1-parameter. 

Here we use Monte Carlo simulations to determine the range 
of potential runaway velocities, orbital separations and merger 
timescales as predicted by the standard model for the formation 
of neutron star/black hole binaries. We include the many recent de- 
velopments in the field which severely affects the input physics and 
the outcome - especially the helium star evolution and the binding 
energy during the common envelope evolution. Also the question 
of kicks imparted to newborn black holes has been revisited after 
the recent detection of the space velocity of the black hole binary 



GRO J1655-40 (Mirabel et al. 2002). Our binaries are launched 
into galaxies with various masses. Initially the binaries are placed 
in starforming regions and we keep track of their motion caused by 
recoil impacts from the supernova explosions of the stellar compo- 
nents. The resulting merger site distributions of the compact bina- 
ries are then explored for the different galactic potentials. 

Population synthesis of massive binaries is also an important 
tool to constrain the local merging rate of NS/BH-binaries in or- 
der to predict the number of events detected by the gravitational 
wave observatories LIGO/VIRGO. As a compact binary continue 
its inspiral, the gravitational waves sweep upward in frequency 
from about 10 Hz to 10 3 Hz, at which point the compact stars 
will collide and coalesce. It is the last 15 minutes of inspiral, with 
~ 16 000 cycles of waveform oscillation, and the final coalescence, 
that LIGO/VIRGO seeks to monitor. LIGO I and LIGO II (the ad- 
vanced interferometer) are expected to detect NSNS inspiral events 
out to a distance of ~ 20 Mpc and ~ 300 Mpc, respectively, accord- 
ing to recent estimates (Thorne 2001). However, as a result of their 
larger chirp masses, merging double black hole binaries (BHBH) 
are detected out to a distance which is ~5 times larger (see e.g. Si- 
pior & Sigurdsson 2002). We have therefore extended our work to 
also include these systems. 

It is the hope that future observations will reveal a combined 
gravitational wave and GRB source. This should be possible given 
the angular resolution of ~ 1 degree on the sky from simultaneous 
detections using both LIGO and VIRGO. To establish such a con- 
nection between these two events would reveal interesting new in- 
sight to the physics involved. Also by comparing the arrival times 
of the gravitational waves and the earliest gamma rays, it should 
be possible to measure the relative propagation speeds of light and 
gravitational waves to an accuracy ~ 10~ 17 (Thorne 2001). 

In our binary notation we distinguish between BHNS, where 
the black hole is formed before the neutron star, and NSBH where 
the black hole is formed after the neutron star. NSNS and BHBH 
refer to double neutron star and double black hole systems, re- 
spectively. In Section 2 we describe our model for making neu- 
tron star/black hole binaries via interacting binary evolution. Our 
modelling of galactic potentials and binary star formation rates are 
presented in Section 3. In Section 4 we summarize our results and 
in Section 5 we follow up with discussions and compare with pre- 
vious work. Finally, our conclusions are given in Section 6. 



2 FORMATION OF NEUTRON STAR/BLACK HOLE 
BINARIES 

In Fig. 1 we have shown the standard evolutionary sequences lead- 
ing to the formation of a double neutron star system. The (recycled) 
neutron stars are detected as radio pulsars. The scenario for pro- 
ducing a black hole/neutron star binary system is similar. As can be 
seen from the figure the progenitor systems evolve through a high- 
mass X-ray binary (HMXB) phase. The formation of a HMXB re- 
quires two relatively massive stars (> 12 M Q ) on the ZAMS. Alter- 
natively, the secondary star can be less massive initially, as long as 
it gains enough material from the primary star so that it will later 
end up above the threshold mass for undergoing a supernova ex- 
plosion (like the primary star). The first mass transfer phase, from 
the primary to the secondary star, is usually assumed to be dynami- 
cally stable (semi-conservative) since the mass ratio, at the onset of 
the RLO, is not too extreme. However later on, when the secondary 
star evolves, all HMXBs must end up in a common envelope phase, 
as the neutron star (or low-mass black hole) is engulfed by the 
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Figure 1. Cartoon depicting the formation of a HMXB and finally a double 
neutron star (and/or black hole binary system). Such a binary will experi- 
ence two supernova explosions. It is always the recycled pulsar which is 
observed in a double pulsar system as a result of its very long spin-down 
timescale compared to a young pulsar (a factor of ~ 10 2 ). If the final system 
is very tight it will coalesce as a result of gravitational wave radiation. 

extended envelope of its companion, in an orbit which is rapidly 
shrinking due to significant loss of orbital angular momentum. This 
stage of binary evolution is extremely important, since the outcome 
is a huge reduction of the orbital separation or merging of the stel- 
lar components - see Sect. 12.41 Stellar winds of massive stars, as 
well as of naked helium cores (Wolf-Rayet stars), are also some of 
the most uncertain aspects of the modelling of HMXB evolution. 
Finally, the physical conditions which determine the formation of a 
neutron star versus a black hole are also still quite unknown. Below 
we shall briefly outline our model of binary evolution. 

2.1 Stellar evolution models 

We used Eggletons numerical stellar evolution code to collect 
a large array of grid points for stars with masses in the inter- 
val 0.9 - 100 M . This code uses a self-adaptive, non-Lagrangian 
mesh-spacing which is a function of local pressure, temperature, 
Lagrangian mass and radius. It treats both convective and semi- 
convective mixing as a diffusion process and finds a simultane- 
ous and implicit solution of both the stellar structure equations 
and the diffusion equations for the chemical composition. New im- 
provements are the inclusion of pressure ionization and Coulomb 
interactions in the equation-of-state, and the incorporation of re- 
cent opacity tables, nuclear reaction rates and neutrino loss rates. 
The most important recent updates of this code are described in 
Pols et al. (1995, 1998) and some are explained in Han et al. 
(1994). In our standard model we used a chemical composition 
of (X=0.70, Z=0.02) and assumed a mixing-length parameter of 
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Figure 2. Evolutionary tracks for a few selected stars in our code (X=0.7(), 
Z=0.02). For the massive stars (12 - 100 M Q ) the final mass, after stellar 
wind mass loss, is written at the end of the tracks. 

a = II H p = 2.0. The uncertainty of these values (and their tempo- 
ral evolution) is not very significant for the results of our work in 
consideration here, given the many major uncertainties of parame- 
ters governing the mass-transfer phases, spiral-in, helium star wind, 
supernova explosion etc. Convective overshooting is taken into ac- 
count in the same way as in Pols et al. (1998) using an overshooting 
constant S m = 0.10. We used de Jager's work (de Jager et al. 1988; 
Nieuwenhuijzen & de Jager 1990) to estimate the wind mass loss 
prior to the mass-transfer phase. Examples of our stellar evolution- 
ary tracks are shown in Fig.|2| We calculated ~ 40000 grid points 
of fine mesh spacing in the HR-diagram in order to perform fast 
Monte Carlo population synthesis. 

2.2 A model of binary evolution 

In general terms our model follows the standard model for the 
formation and evolution of HMXBs. Therefore, in the following 
we shall mainly discuss the issues below where our model ap- 
plies different input physics compared to earlier work. For a gen- 
eral discussion of close binary evolution we refer to, for example, 
van den Heuvel (1994) and Tauris & van den Heuvel (2003). 

To simulate the formation of massive double degenerate bina- 
ries we assume that the initial binary system consists of two zero- 
age main sequence (ZAMS) stars, and chose primary masses in 
the interval 10 - 100 M Q using a Salpeter-like IMF, N(m) oc f?r 2 - 70 
(Scalo 1986). The secondary stars (less massive than their primary 
companions) were chosen in the interval 4 - 100 M Q according to 
the mass-ratio function: f(q) = 2/(1 + qj 1 (Kuiper 1935). We adopt 
the term "primary" to refer to the initially more massive star, re- 
gardless of the effects of mass transfer or loss as the system evolves. 
The initial orbital separations were chosen between 5 - 10000i? Q 
from a distribution flat in log a. 

2.3 Helium stars, SN kicks and remnant masses 

For low-mass stars, the mass of the helium core in post main- 
sequence evolution is practically independent of the presence of an 
extended hydrogen-rich envelope. However, for more massive stars 
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Figure 3. Evolutionary tracks for selected helium stars in our code (Y=0.98, 
Z=0.02). The final mass (after stellar wind mass-loss) is written at the end of 
the tracks. The expansion of these (low-mass) helium stars in close binaries 
often results in a secondary mass-transfer phase. The more massive Wolf- 
Rayet stars do not expand very much during their evolution. Almost all the 
He is burned into C, O and subsequently Ne. After O. Pols (priv. comm.) 

(> 2 M Q ) the evolution of the core of an isolated star differs from 
that of a naked helium star (i.e. a star which has lost its hydrogen 
envelope via mass transfer in a close binary system). It is very im- 
portant to incorporate the giant phases of helium star evolution. Of 
particular interest are the low-mass helium stars (M Hc < 3.5 M Q ) 
since they swell up to large radii during their late evolution (e.g. 
Habets 1986). This may cause an additional phase of mass transfer 
from the naked helium star to its companion (often referred to as 
case BB mass transfer). Recent detailed studies of helium stars in 
close binaries have been performed by Dewi et al. (2002), Dewi & 
Pols (2003) and Ivanova et al. (2003). It should be noticed that he- 
lium cores in binaries have tiny envelopes of hydrogen (< 0.01 M Q ) 
when they detach from Roche-lobe overflow (RLO) mass transfer. 
This tiny envelope seems to have important effects on the subse- 
quent radial evolution of the helium star (e.g. Han et al. 2002). 
The evolution of more massive helium stars (Wolf-Rayet stars) 
is also quite important. There is currently not a clear agreement 
on the rate of intense wind mass-loss from Wolf-Rayet stars (e.g. 
Wellstein & Langer 1999; Nugis & Lamers 2000; Nelemans & 
van den Heuvel 2001). A best estimate fit to the wind mass-loss rate 
of Wolf-Rayet stars is, for example, given by Dewi et al. (2002). We 
have adopted their helium star models (Z=0.02, Y=0.98), including 
wind mass-loss, as calculated by Onno Pols (private communica- 
tion) - see Fig. [3] The uncertainty in determining the wind mass- 
loss rate also affects the threshold mass for core collapse leading 
to a black hole (Schaller et al. 1992; Woosley, Langer & Weaver 
1995; Brown, Lee & Bethe 1999). However, it may well be that 
core mass is not the only important factor to determine the out- 
come. Magnetic field and spin of the collapsing core could also 
play a major role (Ergma & van den Heuvel 1998; Fryer & Heger 
2000). Furthermore, it seems clear from observations that there is 
an overlap in the mass range for making a neutron star versus a 
black hole. In our code we adopted the threshold values given in 
TableQ The threshold mass for the CO-core corresponds to a he- 
lium star which has experienced another RLO (case BB) in its giant 
phase and thereby lost its helium envelope. 



Table 1. Threshold masses for producing black holes and neutron stars 
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10 ~ 12 M * 



* Depending on the evolutionary status at the onset of RLO (case C or B). 
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Figure 4. A momentum kick is assumed to be imparted to all newborn 
neutron stars. The kick magnitude of the asymmetric SN was drawn from 
a weighted sum of the three Maxwellian distributions shown above. The 
direction of the kicks is assumed to be isotropic. Black holes were assumed 
to be accompanied with a smaller kick velocity - see text for description. 



We assumed all neutron stars to be formed in asymmetric SNe 
with a mass of 1.3 M B and chose the 3-D magnitude of the kicks 
from a weighted sum of three Maxwellian distributions of speeds 
with cr„. = 30, 175 and 700 km s" 1 (5%, 80% and 15%, respec- 
tively - see Fig.[4j- The choice of this distribution is partly moti- 
vated by the study of Cordes & Chernoff (1998), but also includes 
a low-velocity component to account for the many neutron stars re- 
tained in rich globular clusters like 47 Tuc (e.g. Pfahl et al. 2002). 
Furthermore, we assumed the direction of the kicks to be isotropic 
(see however Lai, Chemoff & Cordes 2001) and neglected any im- 
pact on the companion star from the ejected shell. The recent de- 
termination of the run-away velocity ( 1 1 2 ± 1 8 km s~ ' ) of the black 
hole binary GRO J 1655-40 (Mirabel et al. 2002) seems to suggest 
that also a stellar collapse leading to the formation of a black hole 
is accompanied by a momentum kick. In our standard model we 
apply a reduced kick in the formation process of black holes. We 
do this by assuming that a similar linear momentum, comparable 
to that imparted to neutron stars, is given to black holes, such that 
their resulting kick velocities are reduced by a factor of M BH /M NS 
with respect to our kick distribution applied on newborn neutron 
stars. For the masses of the black holes we assumed a mass loss 
fraction of 30% during the stellar collapse (Nelemans, Tauris & 
van den Heuvel 1999). 
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2.4 Common envelope and spiral-in evolution 



Except for a few cases in which the two stars in a massive binary 
are born on the ZAMS with almost equal masses (Brown 1995), the 
system will evolve through a HMXB -phase and subsequently en- 
ter a common envelope and spiral-in phase. In our code we follow 
the treatment of common envelope evolution introduced by Dewi 
& Tauris (2000). The same method has been applied recently by 
others, e.g. Podsiadlowski, Rappaport & Han (2002). Hence, we 
calculate the binding energy of the stellar envelope by integrating 
through the outer layers of the donor star: 



Mcore 



GM(r) 



dm + ttti. 



Mo 



Udm 



(1) 



where the first term is the gravitational binding energy and U is the 
internal thermodynamic energy (Han et al. 1994, 1995). The value 
of a th depends on the details of the ejection process, which are very 
uncertain. A value of equal to or 1 corresponds to maximum 
and minimum envelope binding energy, respectively. By adopting 
the energy formalism of Webbink (1984) and de Kool (1990), we 
then equate this envelope binding energy with the release of orbital 
energy in order to estimate the reduction of the orbit during spiral- 
in: £bind = rjcE A£ or b where tjce describes the efficiency of ejecting 
the envelope, i.e. of converting orbital energy into the kinetic en- 
ergy that provides the outward motion of the envelope. The total 
change in orbital energy is then simply given by: 



A£„, 



GM corc Mi GM ionot Mi 



(2) 



2 Of 2 a\ 

where M colc = Md onor _ M cm is the mass of the helium core of the 
evolved donor star; Mi is the mass of the companion star; a, is the 
initial separation at the onset of the RLO and a { is the final orbital 
separation after the CE-phase. The orbital separation of the surviv- 
ing binaries is quite often reduced by a factor of ~ 100 as a result 
of the spiral-in. If there is not enough orbital energy available to 
eject the envelope the stellar components will coalesce. For mas- 
sive stars (M > 10 M Q ) this method results in so-called /1-values 
< 0.1 - 0.01 (Dewi & Tauris 2001; Podsiadlowski, Rappaport & 
Han 2002). These values are much smaller than the constant value 
of A ~ 0.5 often used in the literature (e.g. Hurley, Tout & Pols 
2002; Sigurdsson & Pols 1999; Fryer, Woosley & Hartmann 1999; 
Portegies Zwart & Yungelson 1998, or alternatively, A rj C E is used as 
a single free parameter: Belczynski, Kalogera & Bulik 2002; Bel- 
czynski, Bulik & Kalogera 2002; Perna & Belczynski 2002). This 
result has the very important consequence that in our study many 
HMXBs will not produce double neutron star/black hole systems 
because they merge during their subsequent CE-phase (however, 
those systems that do survive are generally in tight orbits often 
leading to additional helium star mass transfer and formation of 
double degenerate systems merging on a short timescale). 
It should be noticed that the exact determination of A depends on 
how the core boundary is defined (see Tauris & Dewi 2001 for a 
discussion). For example, if the core boundary (bifurcation point 
of envelope ejection in a CE) of the 20 M e in Fig|5|is moved out by 
0.1 M Q then A is typically increased by a factor of ~2. 

2.4.1 Helium star common envelope 

We take into account that some helium stars, in their giant stages, 
regain contact with their compact companion star. This leads to a 
second mass-transfer phase resulting in a naked CO-core, if the 
stellar components do not coalesce. From detailed helium star anal- 
ysis, Dewi & Pols (2003) discovered that there is a critical mass 
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Figure 5. The /1-parameter for a 20 M Q star as a function of stellar radius. 
The upper curve includes internal thermodynamic energy (a t h = 1 ) whereas 
the lower curve is based on the sole gravitational binding energy (a^ = 
0) - see Eq. [3 There is a factor ~2 in difference between the /t-curves in 
accordance with the virial theorem. It is a common misconception to use a 
constant (and large) value of A =s 0.5 marked by the dashed line. In this plot 
the stellar core is defined as the region where X< 0.10. See text. 



(M™ 1 ~ 3.3 M Q ) above which the RLO from the helium star to 
a neutron star is dynamically stable. Below this value the system 
will evolve into a CE and spiral-in phase as a consequence of the 
rapid expansion of the low-mass helium star. We adopted this value 
in our code to determine the fate of helium star CE (see however 
Ivanova et al. 2003 for a different point of view). Typical helium 
star /1-values are 0.01 - 0.2 depending on radius and a^. 



2.4.2 Released accretion energy during CE evolution 

Even though the spiral-in timescale is very short (< 1000 yr) 
for common envelope evolution (Taam & Sandquist 2000), and 
the accretion rate is limited by the photon radiation pressure (the 
Eddington limit), the gravitational potential energy release from 
the ~ 10~ 5 M accreted onto the neutron star or black hole con- 
tribute to expel the envelope of the donor star. We modelled this 
effect by assuming a total energy release of A£ acc ~ T ce L Edd = 
3.8 x 10 48 (Mx/M ) erg, and twice this amount for helium star CE, 
where M\ is the mass of the compact object. Given the high col- 
umn density inside the common envelope a significant part of this 
energy will be absorbed and facilitate the ejection of the envelope 
and thereby increase the post-CE orbital separation by a factor: 



at (1 + A£ acc /£ bmd ) 



(3) 



The average value of this factor is ~ 1.8 for all our systems evolv- 
ing through a CE-phase (the typical value is only ~ 1.3, but a few 
systems with very evolved donors have A£ acc > £bind and a corre- 
sponding correction factor up to about 5). 



2.5 RLO and the orbital angular balance equation 

The orbital angular momentum of a binary system is given by: 



\fXp\ : 



M 



■Sla 2 VT 



(4) 
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where a is the separation between the stellar components, M\ and 
M 2 are the masses of the accretor and donor star, respectively, 
M = Mi + M 2 and the orbital angular velocity, SI = ^GM/a 3 . 
Here G is the constant of gravity. Tidal effects acting on a near- 
RLO (giant) star will circularize the orbit on a short timescale of 
~ 10 4 yr (Verbunt & Phinney 1995) and in the equation below 
we therefore disregard any small eccentricity (e = 0). A simple 
logarithmic differentiation of the above equation yields the rate of 
change in orbital separation: 



i=2 — -2 — -2^- + M[ +Ml ~ 
a J orb Mi M 2 M 

where the total change in orbital angular momentum is: 



J orb 
"orb 



^ -Anb ^ *As ^ "ml 
J orb "orb "orb J orb 



(5) 



(6) 



These two equations constitute the orbital angular momentum bal- 
ance equation. The first term on the right-hand side of this equation 
gives the change in orbital angular momentum due to gravitational 
wave radiation (Landau & Lifshitz 1958). The second term in Eq.[6] 
arises from so-called magnetic braking (this mechanism is not rel- 
evant for HMXBs). The third term (ji s /J OI b) on the right-hand side 
of Eq. {6j describes possible exchange of angular momentum be- 
tween the orbit and the donor star due to its expansion or contrac- 
tion (Tauris 2001). Finally, the last term in Eq. |6j represents the 
change in orbital angular momentum caused by mass loss from the 
binary system. This is usually the dominant term in the orbital an- 
gular momentum balance equation and its total effect is given by: 



Jmi _ a+Pq 2 +8y(\ + qf M 2 
J or b l+q M 2 



(7) 



where q = M 2 jM\ is the mass ratio, and where a, p and r5 are 
the fractions of mass lost from the donor in the form of a di- 
rect fast wind, mass ejected from the vicinity of the accretor and 
from a circumbinary coplanar toroid (with radius, a r = y 2 a), re- 
spectively - see van den Heuvel (1994) and Soberman, Phinney & 
van den Heuvel (1997). The accretion efficiency of the accreting 
star is thus given by: e = 1 - a — p — 6, or equivalently: 



dMi = -(1 -a-p-6)dM 2 



(8) 



where dM 2 < (Mi refers to the donor star). These factors will 
be functions of time as the binary system evolves during the mass- 
transfer phase. The general solution for calculating the change in 
orbital separation during the X-ray phase is found by integration of 
the orbital angular momentum balance equation (Eq.[5}: 



qo 



-a-p+5 



q+l\ 




qo + ij 





(9) 



ac 2 +J3+ r 6(\-€)- 

eq + 1 \ 



eq + 1 



where the subscript '0' denotes initial values and T;, is factor of 
order unity to account for tidal spin-orbit couplings ('i s ). 

fn our standard model for the RLO we assumed a direct wind 
mass-loss fraction of a = 0.20 and neglected the possibility of any 
circumbinary toroid (i.e. 6 = 0). Furthermore, we assumed the ac- 
cretion rate onto the compact object to be limited by the Edding- 
ton luminosity and hence: p = max((\M Aovim \ - M Ec | d )/|M donor |, o), 
where |M donor | is the mass-transfer rate from the donor star. For a 
typical neutron star accreting hydrogen MecM - 1.5 x 10~ 8 M e yr -1 . 



2.6 Gravitational waves and merging NS/BH binaries 

In the Newtonian/quadrupole approximation (a <k A gwr ) the merg- 
ing timescale of a compact binary with semi-major axis, a and ec- 
centricity, e is given by Peters (1964): 



T swr (a ,e ) ■ 



1? C* H 
19 p Jo 



«0 „29/19 



[1 +(121/304)e 2 ] 



2-1 1 181/2299 



(1 _ e 2)3/2 



-de (10) 



where 



a(l-e 2 ) 

,,12/19 



[1 +(121/304)e 2 ] 



,2-1-870/2299 



(ID 



can be found from the initial values (a ,eo)> an d the constant p is 
given by: 

„ 64G 3 



5c 5 



(12) 



Here M and n denote the total mass and the reduced mass of the 
system, respectively. The integral given above must be evaluated 
numerically. However, for circular orbits the merging timescale can 
easily be found analytically: t"™. = a* /Ap. The timescale is very 
dependent on both a and e. Tight and/or eccentric orbits spiral-in 
much faster than wider and more circular orbits - see Fig.[l0l For 
example, we find that the double neutron star system PSR 1913+16 
(P orb = 7.75 hr, M NS = 1.441 and 1.387 M Q , respectively) with an 
eccentricity of 0.617 will merge in 302 Myr; if its orbit was circular 
the merger time would be five times longer: 1.65 Gyr! 



3 STAR FORMATION RATES, GALACTIC MODELS 
AND STARFORMING REGIONS 

We have assumed that the star formation rate has been con- 
tinuous in the Galactic disk as inferred from observations (e.g. 
Gilmore 2001). Furthermore, following Hurley, Tout & Pols (2002) 
we assume that one binary (Mi > 0.8 M Q ) is born in the Galaxy per 
year. Hence, for the formation rate of a massive binary with primary 
star mass > 10 M Q and secondary star mass > 4 M Q we simply use: 



$bsfr = /0.01 yr 



(13) 



where / is a factor ~ 1 to scale all our derived formation and merger 
rates. We fix this rate over the lifetime of the Galaxy (~ 12 Gyr). 
This prescription should facilitate comparison with other studies. 



3.1 Galactic potentials 

We have evolved our binaries in the potential of spiral galaxies. 
The adopted potential is a model for the Milky Way made by Flynn, 
Sommer-Larsen & Christensen (1996) with a few corrections to the 
parameters (Sommer-Larsen, private communication 2002). The 
potential <t>(R, z) is modelled in cylindrical coordinates: R being the 
distance in the plane from the center of the galaxy, and z being the 
height above the plane. The potential is a sum of a dark matter halo 
0#, a central component <t> c and a disk O d . The dark matter halo 
is assumed to be spherically symmetric: 

f„ = l -V 2 H \n(r 2 + rl) (14) 

where r 2 = R 2 + z 2 ■ The central part of the potential is made-up of 
two spherical components: 

GM Cl GM Cl 



(15) 
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Table 2. Adopted parameters for the Galactic potential 



Component 

Bulge/Stellar halo 
Central Component 
Disk 




BHBH BHNS NSBH NSNS 



Here G is the gravitational constant, Mc, and f*c, are the mass and 
the core radius of the bulge/stellar halo term, while Mc 2 and rc 2 
are the mass and core radius of the inner core. The disk potential is 
modelled using a combination of three Miyamoto-Nagai potentials, 
see Miyamoto & Nagai (1975): 

-GM D 

= — " . n= 1,2,3 (16) 

^(R 2 + [a„ + VPTFp) 

The parameters b and a„ are related respectively to the scaleheight 
and scalelength of the disk, and M Dn are the masses of the three disk 
components. The parameter values for the Milky Way are found in 
Tabled 

To obtain potentials for galaxies with other masses than the 
Milky Way, we have scaled the parameters with the square root of 
the galaxy mass (assuming that the ratio of the components of the 
galaxies remains constant). This gives a constant surface brightness 
as observed, see e.g. Binney & Tremaine (1994). 

Star formation takes place in two very distinct physical en- 
vironments. Mostly in the extended disks of spiral and irregular 
galaxies, but also a significant fraction in the compact dense gas 
disks in the centers of galaxies. The central star formation happens 
in starbursts with very high star formation rates for 0.1-1.0 Gyr 
(Kennicut 1998), while star formation in the extended disk is steady 
or slowly exponentially decreasing (Bruzual & Chariot 1993). We 
have assumed that the distribution of birth places follow the den- 
sity distribution of the extended disk. The latter can be found by 
applying Poisson's equation to the disk potential: 

V 2 = AnGp D (17) 

Finally, we assume that the binaries are born with a velocity equal 
to the local rotational velocity. 

4 RESULTS 

The results presented in this paper are based on a typical simulation 
of 20 million ZAMS binaries with standard parameters as given in 
Tabled Of these systems 99.9 % did not make it to the final stage. 
Most of these systems either merged during the CE and spiral-in 
phase or were disrupted by a supernova explosion. 



Figure 6. Relative formation ratios of binaries with neutron stars and/or 
black holes which will coalesce within 10 Gyr. 

4.1 Relative formation ratios and merging timescales 

The aim of this paper is to focus on GRB progenitors and burst 
sources of gravitational wave radiation. Therefore, we start out by 
showing in Fig. [6| the relative formation ratios of BHBH, BHNS, 
NSBH and NSNS systems which will merge within 10 Gyr. A total 
of 23 413 such sources is formed over a time interval of ~2 Gyr 
from 20 million massive ZAMS binaries (see Table[3]for input pa- 
rameters). It is seen that double black hole systems (BHBH) dom- 
inate this population and that only few close-orbit NSBH systems 
are formed. However, from Fig.Qit is seen that the relative forma- 
tion ratios of mergers do not represent the overall relative forma- 
tion ratios of all massive double degenerate systems. Many NSBH, 
BHNS and BHBH systems are formed in very wide orbit systems 
(^orb ^ 100 - 1000 yr) which will never merge (r mcrgc ~ 10 21 yr). 
For these three categories the formation of such wide-orbit systems 
is even seen to be dominating. Almost all of the binaries formed in 
this right-hand peak of the bimodal T metge distribution are caused by 
binaries which did not experience RLO from their secondary star 
(causing the orbits to widen subsequently from stellar wind mass 
loss and the effects of the second SN explosion). The majority of 
the wide-orbit BHNS systems also evolved without primary star 
RLO and thus remained detached binaries during their entire evolu- 
tion. The same condition applies to ~ 1/3 of the wide-orbit BHBH 
systems, but none of the NSBH systems. The reason for the latter is 
that without RLO the secondary star will never end up with a mass 
above the threshold limit for making a black hole. Another interest- 
ing feature is that almost half of all the BHNS systems which merge 
within 10 Gyr evolve without primary star RLO. However, in these 
systems the black hole is shot into a closer orbit (due to the kick) 
which later, as a result of nuclear expansion, causes the secondary 
star to fill its Roche-lobe. This gives rise to the formation of a com- 
mon envelope and spiral-in phase leading to a close final orbit. The 
same scenario is also responsible for ~ 5 % of the close-orbit NSNS 
systems. This formation channel (direct-supernova) of close-orbit 
systems was first discussed by Kalogera (1998) in the context of 
LMXBs. Another comment to Fig.0is that the population of close- 
orbit NSBH systems is clearly seen to be subdivided into two dis- 
tinct populations. The systems with T mcrec < 9 Gyr evolved through 
helium star RLO (case BB) onto the NS whereas the systems with 
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Figure 7. Histograms of merging timescales for binaries with neutron stars and/or black holes. The populations of NSBH, BHNS and BHBH are seen to be 
bimodal with a large component of very wide binaries with T merge ~ 10 21 yr. See text for a detailed discussion. 




Figure 8. Histograms of merging timescales for binaries with neutron stars and/or black holes which coalesce within 1 Gyr. 
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Table 3. Adopted standard parameters for our binary population synthesis code. A typical simulation is based on 20 million ZAMS binaries. 



Parameter 


Symbol 


Value 


Note 


JVLass 01 primary star 




10—1 on M 


Snlnptpr IMP ( 9 1\ 




M 


4 - M 


thp Hi striHiitinn is Hrnwn to match thp n-ntctnhntinfl 

UIL LLI>IJ JL'LILIV'II 1> U1LIW11 HI lllCtLL.ll L11L* LI LI J .M 1 J U 11 1 H 1 1 J 


Initial ZAMS separation 


ao 


5- IOOOORq 


a flat distribution in log a 


Stellar wind mass loss 


a 


0.20 




Timescale for RLO 


Afrlo 


~ ^ ^"thermal 


generally a good estimate 


Critical mass ratio for CE 


<Zee 


2.5 


for a hydrogen donor star (at onset of RLO after wind mass loss) 


Timescale for CE 


Afce 


1000 yr 


upper limit 


Efficiency of CE 


Ice 


0.5 


+ liberated ii acc , see Eq.l3l 


A for CE 


A 


0.006 ~ 0.4 


individually calculated from stellar structure calculations, see Sec. 12. 41 


Core boundary parameter 


fx = r/A 


2 


since X ^ 0.10 may often underestimate the core mass (Tauris & Dewi 2001) 


Critical He-star mass for CE 


<' 


3.3 Mq 


from recent helium donor star calculations (Dewi & Pols 2003) 


,1 for He-star CE 


Aae 


0.1 


typical value from (Dewi & Pols 2003) 


Maximum accretion rate 


M max 


M E dd 


an Eddington accretion limit is applied to both NS and BH 


kick-distribution for NS 


w 


- 2000 km s" 1 


3 comp. Maxwellian [cr=30 (5%), 175 (80%) and 700 km s _1 (15%), respectively] 


kick-distribution for BH 


»'bh 


- 2000 km s _1 


same as above, but scaled: = w(M ni /M^) to yield same linear momentum 



Table 4. The formation rates, the merger rates and the present population of 
massive double degenerate binaries in the Galactic disk - see text. 



Systems 


Formation rate 


Merger rate 


/VGal 


NSNS 


1.8 x 10~ 6 yr-' 


1.5 X 10~ 6 yr 1 


3 700 


NSBH 


6.3 x ur 7 yr 1 


8.4 X 10~ 8 yr 1 


6400 


BHNS 


l.i x ur 6 yr 1 


5.0 X 10~ 7 yr 1 


6900 


BHBH 


3.7 X I0~ 5 yr" 1 


9.7 X lO^yr -1 


320000 



9 < T mcrge /Gyr < 12 evolved without helium star mass-transfer. 
A close-up of the distribution of merging timescales < 1 Gyr is 
presented in Fig. [8] The various populations have rather different 
distributions of r mcrg e - ranging from extremely short (NSBH) to 
somewhat fiat (BHNS). 

4.2 Galactic formation rates, merger rates and the present 
population of massive double degenerate binaries 

In Table|4]we list our best estimate for the Galactic formation rate, 
the merger rate and the present population of massive double de- 
generate binaries. The present population in the Galactic disk is a 
rough estimate based on the difference in formation rate and merger 
rate integrated over the age of the Milky Way (12 Gyr) and cor- 
rected for: I) the systems escaping the gravitational potential of our 
Galaxy (< 10%) and II) the contribution from the current popu- 
lation of potential mergers recently formed. The formation rate of 
NSNS systems is larger than that of both NSBH and BHNS sys- 
tems. However, as we noticed from Fig.0 discussed in the previ- 
ous section, these latter systems form a large fraction of very wide- 
orbit binaries which will never merge and therefore accumulate in 
numbers with the age of the Galaxy. This is also seen from their rel- 
atively small merger rates. It is thus expected that the Galaxy hosts 
a larger number of NSBH and BHNS systems compared to NSNS 
systems. The BHBH systems, however, are seen to dominate the 
overall Galactic population of massive compact binaries. We es- 
timate a present population above 300 000 systems in the Milky 
Way. Although such systems cannot be detected observationally yet 
(since they do not emit electromagnetic radiation) we may hope to 
detect continuous gravitational waves from some of these systems 
in our backyard once LISA has been launched. 



In order to understand the large formation rate of BHBH bi- 
naries, and the origin of the different massive compact binaries, it 
is usefull to look at Fig. [5] It is clearly seen that the initial param- 
eter space for the stellar component masses favour the formation 
of BHBH binaries. The lines separating the different populations in 
the diagram can largely be understood from the threshold masses 
listed in Table Q The dividing lines are "fuzzy" since the initial 
orbital separation is also a determining factor for the outcome of 
the binary stellar evolution (i.e. a helium star's mass depends on at 
what evolutionary stage its progenitor star lost its hydrogen enve- 
lope). The lower curved line in the plot marks the onset of common 
envelope evolution. The line is curved (rather than straight) since 
the mass-ratio limit q,. c is evaluated after the wind mass loss is ac- 
counted for, prior to the RLO, and more massive stars loose a rela- 
tively larger fraction of their mass in the form of a stellar wind (see 
Fig-12}- We notice that no BHBH systems are formed which evolved 
through a CE-phase during the mass transfer from the primary star. 
The reason for this is the small values of the /1-parameter obtained 
from real stellar structure models. This reflects that the envelopes 
of these massive stars are tightly bound to their cores so that the bi- 
naries do not have enough orbital energy to expel the envelopes in 
the spiral-in phase. Therefore the systems will merge (and perhaps 
form Thorne-Zytkow objects). 

4.3 Galactic double neutron star (NSNS) binaries 

Simulating the formation of double neutron star binaries is a cru- 
cial test for the computer code and the physics assumptions be- 
ing used. In Fig. 1101 we have plotted the distribution of NSNS 
binaries at birth in a (separation, eccentricity)-diagram. We have 
partly used the location of the four known Galactic NSNS sys- 
tems in the disk (see Table|5J to constrain our input parameters so 
we can reproduce these observed systems from simulations. Espe- 
cially PSR J15 18+4904 (Nice, Sayer & Taylor 1996) is a challenge 
for any binary population synthesis code given its low eccentric- 
ity (0.27) and relatively wide orbit P or \, = 8.6 days {a ^ 25 R ). 
We find that our simulated population fits the observations well 
given that the majority of the simulated NSNS systems merge on 
a timescale much shorter than the lifetime of the observed (recy- 
cled) pulsars. Furthermore, NSNS systems with a < 1 R Q (P„rb £: 
100 min.) are very difficult to detect as a result of the large acceler- 
ation of their pulsed signals. It is possible that this plot may help to 
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Figure 10. The simulated distribution of newborn double neutron star binaries in the (separation, eccentricity)-plane. Isochrones for the merging time of the 
double neutron star binaries are also shown. We calculated these assuming two 1.4 M Q neutron stars in each binary. The curves correspond to values of: 
3 X 10 5 yr, 3 Myr, 30 Myr, 300 Myr and 3 Gyr (from left to right), respectively. The four NSNS systems detected in the Galactic disk are indicated with stars. 



Table 5. The four observed double neutron star systems in the Galactic disk 
(for references see Nice, Sayer & Taylor 1996 and Lyne et al. 2000). 



NSBH 




BHNS 

i i i i i 

60 80 



100 



M p /M 



Figure 9. Initial primary and secondary stellar masses for 20000 ZAMS 
binaries evolving all the way to form massive double degenerate systems. 



converge some of the many various input parameters used in popu- 
lation synthesis codes today. 

It is instructive to look at the distribution of kicks imparted 
to neutron stars in surviving binaries which evolved all the way 
to the final NSNS stage. This is shown in Fig. 1111 The difference 
between the kicks imparted during the first and the second SN ex- 
plosion is very significant: the average values of the selected kicks 
are W[ = 53 km s and w% = 263 km s , respectively. In order 



PSR-name 


sep. 


ecc. 


M = M\ + M2 


T = P/2P 




Re 




M Q 


Myr 


J 1518+4904 


24.9 


0.249 


2.69 


20000 


B 1534+12 


3.33 


0.274 


2.68 


200 


J 1811-1736 


41.9 


0.828 


2.6 


900 


B 1913+16 


2.81 


0.617 


2.83 


100 


B 2127+1 1C* 


2.86 


0.681 


2.71 


100 



This system resides inside a globular cluster and has a different origin. 



for the systems to avoid merging in the CE and spiral-in phase of 
the evolved secondary stars, the systems must be in a fairly wide 
orbit prior to the HMXB and subsequent CE and spiral-in phase. 
This means that these systems can only survive very small kicks in 
the first SN explosion in order to avoid disruption. On the other 
hand, after the CE and spiral-in phase the orbits of the surviv- 
ing systems will always be tight and therefore the systems will be 
able to withstand even very large kicks in the second SN explo- 
sion. Actually, the distribution of selected kicks during the second 
SN resembles the input trial distribution of kicks plotted in Fig. [5] 
The corresponding systemic recoil velocities after the SNe were 
vi = 8.6 km s _I and v 2 = 220 km s -1 , respectively. Thus, when 
considering the runaway properties of NSNS systems it is sufficient 
to consider the effects from the second SN. 

To get an idea of the distance travelled by NSNS binaries from 
their birth place prior to their coalescence, one simply has to mul- 
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Figure 12. Distribution of newborn double neutron star binaries with respect to orbital separation, merging timescale, systemic (recoil) velocity from the 
second SN explosion and the maximum distance from the origin of birth in a host galaxy - see text for details. 
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Figure 11. The distribution of kicks imparted to newborn neutron stars in 
sun'iving binaries. The difference in selected kicks between the first and 
second SN are expected from a binary evolutionary point of view - see text. 



tiply their space velocity with the merging timescale. The space 
velocity is caused by the recoil of shell ejection combined with 
the asymmetric kick given to the neutron star at birth. In Fie. 1121 
we have plotted histograms of final orbital separation, merger time, 
space velocity and the estimated distance travelled before the merg- 
ing event. In order for other readers to compare with our results we 
have in this plot simplified both the space velocity and the distance 
travelled. For the space velocity shown here it is assumed that the 



binary was at rest prior to the second SN (which is a good approx- 
imation, see above). The other approximation applied here is that 
the binary motion was without influence from the Galactic potential 
(this approximation gets exceedingly worse with increasing merger 
times). Hence, the estimated distance is an upper limit. We see that 
typical systems travel ~10 kpc prior to their merging and that the 
distribution of distances travelled is rather broad. It should be men- 
tioned that in our simulations presented throughout this paper the 
space velocities resulting from both SNe are added as vectors and 
the binary is moving in a Galactic potential as described in Sec- 
tion 3. 



4.4 Offset distribution of mergers from their host galaxies 

Recent work by Bloom, Kulkarni & Djorgovski (2002) shows that 
the median projected angular offset of 20 long-duration GRBs is 
only ~1. 3 kpc (0."17 for an assumed cosmology of fl A =0.7,fl M = 
0.3 and H = 65 km s _I Mpc~'). In Fig. ll3l we have plotted our cu- 
mulative distribution of calculated projected radial distances in var- 
ious host galaxies. For a random orientation of the radial distance 
with respect to the line-of-sight one finds the average projected ra- 
dial distance (in the plane of the sky) by multiplying with a fac- 
tor n/4. The panels in the left column show the result of launching 
our binaries into four galaxies with different masses (0.001-1 times 
the mass of the Milky Way). We see that roughly half of all mas- 
sive compact binaries merge within a projected radial distance of 
4 kpc for a 0. 1 M MW galaxy. Hence, the majority of these systems 
merge inside their host galaxies. However, the median distance is 
larger than the median value observed for the long-duration bursts. 
At a projected radial offset of 1.3 kpc we find that less than 25% 
of all BHNS, NSBH and NSNS systems merge (independent of 
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Figure 13. The panels in the left column shows the cumulative distribution of systems merging within the projected radial galactic distance given on the x-axis. 
Only binaries merging within 10 Gyr are included. We launched our binaries into four galaxies with different masses (0.001-1 times the mass of the Milky 
Way). The panels on the right-hand side are for a 0.1 Mmw galaxy. The curves represent four different merging timescale intervals - see text. 
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Table 6. The expected LIGO/VIRGO detection rates of compact mergers. 




Galaxy age (Gyr) 

Figure 14. Merger rates (relative to present merger rates) as a function of 
the age of our Galaxy for NSNS systems (solid line), BHNS+NSBH sys- 
tems (dashed line) and BHBH systems (dotted line) for a constant Galactic 
star formation rate. The curves show the emergence of a "steady state". 



the mass of the host galaxy). In general, it therefore seems clear 
that merging compact objects most likely cannot be the progenitors 
of long-duration GRBs (as also concluded by Bloom, Kulkarni & 
Djorgovski 2002). However, from a pure kinematical point of view 
one cannot rule out that some of the long GRBs with afterglows 
could be associated with merging compact objects. It is also an in- 
teresting question to what extent such a merging event outside a 
galaxy would produce an afterglow at all. It should be noticed that 
the merging BHBH systems are only shown for comparison - they 
cannot be progenitors of GRBs. Note, that the central panel shows 
the mixed population of both NSBH and BHNS systems. The pan- 
els on in the right-hand side of the figure simply show how the 
cumulative fraction of systems merging increases as a function of 
projected radial galactic distance for a 0.1 M MW galaxy. The curves 
represent four different merging timescale intervals: 0-100 Myr, 
100-500 Myr, 0.5-2.5 Gyr and 2.5-10 Gyr, respectively. We find 
that ~ 10 % of all NSNS systems formed escape our Galaxy. 

Fig. ll4l shows the merger rates of NSNS systems (solid line), 
BHNS+NSBH systems (dashed line) and BHBH systems (dotted 
line) for a constant Galactic star formation rate. Each of the merger 
rates are scaled relative to its present merger rate. In calculating 
this rate, only binaries merging within the age of the Milky Way 
are included. It is seen that the Milky Way today is in a "steady 
state" (i.e. the formation rate of systems in close orbits merging 
within 10 Gyr is equal to the number of systems merging per unit 
time). 

4.5 Prospects for LIGO/VIRGO detection rates 

LIGO I and LIGO II are expected to detect NSNS inspiral events 
out to a distance of ~ 20 Mpc and ~ 300 Mpc, respectively, accord- 
ing to recent estimates (Thorne 2001). This corresponds to wave 
amplitudes of roughly 10~ 20 > h > 10~ 22 . As a result of the much 
larger chirp mass (M chirp = fi 3/5 M 2/5 ) of the BHBH mergers, such 
binaries will be detected out to a distance luminosity, d L oc M 5 l 6 

J ' ^ chirp 

(Finn 1998) which is about 4 times larger in average according to 
our estimates. Hence, the ratio of detected event rates for BHBH 
mergers relative to NSNS mergers is 4 3 ~ 64 times larger than the 



Systems 


Galactic merger rate 


LIGO I 


LIGO II 


NSNS 


1.5 X l(r 6 yr-' 


6.0 X 10~ 4 yr 1 


2.0 yr 1 


NSBH 


8.4 x 1(T 8 yr" 1 


1.7 X lO^yr 1 


0.6 yr 1 


BHNS 


5.0 X ur 7 yr 1 


1.0 X 10~ 3 yr 1 


3.4 yr 1 


BHBH 


9.7 x 10~ 6 yr 1 


2.5 x 10"' yr 1 


840 yr 1 



corresponding ratio of such mergers in our Galaxy. Thus BHBH 
mergers will be completely dominant for LIGO detections (as noted 
by Sipior & Sigurdsson 2002). The exact value of the chirp mass 
ratio depends mainly on the black hole masses estimated at birth 
as well as the possibility of hypercritical accretion onto neutron 
stars. BHNS and NSBH mergers have a chirp mass which is roughy 
2 times larger than that of the NSNS systems and thus the corre- 
sponding detected ratio of such mergers, relative to NSNS mergers, 
is ~5 times larger than the Galactic ratio. 

In order to extrapolate the Galactic coalescence rate out to 
the volume of the universe accessible to LIGO, one can either use 
a method based on star formation rates, galaxy number density 
or a scaling based on the B-band luminosities of galaxies. Using 
the latter method Kalogera et al. (2001) found a scaling factor of 
(1.0 - 1.5) x 10~ 2 Mpc" 3 , or equivalently, ~ 400 for LIGO I (out to 
20 Mpc for NSNS mergers). Since LIGO II is expected to look out 
to a distance of 300 Mpc (NSNS mergers), the volume covered by 
LIGO II is larger by a factor of (300/20) 3 and thus the scaling fac- 
tor in this case, relative to the coalescence rates in the Milky Way, 
is about 1.3 x 10 6 . Therefore, the expected rate of detections from 
galactic field NSNS inspiral events is roughly 2 yr~' for LIGO II 
(see Table|S| and an impressive 840 yr~' for the BHBH mergers! 
Even LIGO I may, with a bit of luck, detect an event from a BHBH 
collision. Beware, that our estimated merger rates should be con- 
sidered as lower limits given that compact mergers in globular clus- 
ters probably also contribute significantly to the total merger rates 
(Portegies Zwart & McMillan 2000). The cosmological implica- 
tions of gravitational wave observations of binary inspiral are also 
interesting to note (e.g. Schutz 1986; Finn 1997). 



5 DISCUSSION 

5.1 The dependence on input parameters 

Monte Carlo simulations is a powerful tool for investigating com- 
plicated interactions which depend on many parameters. Often the 
difficult task is not to produce a reasonable code but to analyse 
the results, pinpoint the important physical parameters behind the 
major trends in a simulated population and try to understand these 
results in terms of relatively simple physics. 

In Table we summarize the effects of various input param- 
eters on the formation and evolution of massive double degenerate 
systems. Model A is our standard model. In all other models we 
modified one parameter at a time, keeping all the other parameters 
as in our standard model. Model B having a constant A = 0.5 is a 
very bad approximation (see Fig. [5}, but nevertheless it is appar- 
ently still used in many computer codes. For example it is clear 
that the formation channel recently presented in Fig. 1 of Bel- 
czynski, Bulik & Kalogera (2002) will not work if one applies a 
reasonable realistic value for the binding energy of the envelope. 
At stage V in their figure, one can use a stellar structure model 
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Table 7. Galactic relative formation ratios, formation rates and merger rates 
of massive compact binaries. Model A is our best estimate. See text. 



Model/ Rel. formation Formation rate Merger rate 

Systems % Myr~' Myr~' 



A) standard 

NSNS 4.4 1.8 1.5 

NSBH 1.5 0.63 0.08 

BHNS 2.7 1.1 0.50 

BHBH 91.3 37 9.7 



B) A = 0.5 constant 

NSNS 7.6 20 17 

NSBH 2.3 6.1 3.7 

BHNS 1.9 4.9 1.3 

BHBH 88.2 230 76 



C) kick: 0.5 oy 

NSNS 6.5 7.2 3.3 

NSBH 3.6 4.0 0.34 

BHNS 3.4 3.8 1.3 

BHBH 86.5 96 18 



D) core: /j = 1 

NSNS 3.6 1.1 0.95 

NSBH 1.6 0.48 0.04 

BHNS 2.9 0.61 0.16 

BHBH 92.7 27 1.3 



E) IMF: / oc M ; 7 2 

NSNS 4.4 1.9 1.4 

NSBH 1.6 0.66 0.08 

BHNS 2.9 1.2 0.54 

BHBH 91.0 37 9.8 



F) wind: a = 

NSNS 2.8 1.1 0.83 

NSBH 3.5 1.4 0.22 

BHNS 2.0 0.77 0.22 

BHBH 91.8 36 7.2 



G) g ce = 2.0 

NSNS 3.6 1.3 1.0 

NSBH 1.5 0.54 0.08 

BHNS 3.7 1.3 0.44 

BHBH 91.2 32 6.9 



H) gee = 3.0 

NSNS 5.2 2.4 1.9 

NSBH 1.5 0.71 0.11 

BHNS 3.7 1.7 1.1 

BHBH 89.5 41 12 



I) ffce.He = 2.5 

NSNS 2.3 0.93 0.54 

NSBH 1.5 0.62 0.02 

BHNS 2.9 1.2 0.57 

BHBH 93.3 38 9.7 



J) A£ acc = 

NSNS 4.2 1.4 1.3 

NSBH 2.2 0.71 0.06 

BHNS 2.2 0.71 0.21 

BHBH 91.4 30 3.4 



average (D— >J) 

NSNS 3.7 1.5 1.1 

NSBH 1.9 0.74 0.09 

BHNS 2.9 1.0 0.47 

BHBH 91.6 34 7.2 



and estimate a value of A = 0.23 for the R = 82 donor star 
of mass 14.1 M Q . This results in a post-CE separation which is 
roughly 3 times smaller than the radius of the descendant helium 
star. Hence, this system is doomed to coalesce (if one excludes the 
internal thermodynamic energy when estimating the value of A the 
problem even exacerbates). Model C seems to give too small val- 
ues of the kick velocity compared to observations of radio pulsars 
(Cordes & Chernoff 1998). At the bottom of the table the average 
result of models D-J is given. It is interesting to notice that these 
values are not very different from our standard model results. 



5.2 The Galactic supernova rate 

We have based the calibration of our models on the star formation 
rate for the Milky Way (see Eq. ll3i . In order to check the validity 
of this calibration we have compared our simulated Galactic super- 
nova rate of type Ib/c SNe with the empirical estimates of Cappel- 
laro et al. (1999). If we assume that all type Ib/c SNe originate from 
the collapse of naked stars which have lost their hydrogen/helium 
envelopes as a result of mass transfer in a binary system, then we 
find a rate of: /0.28 (100 yr)~'. We can compare this result with 
a rate of 0.14 SNu (1 SNu = 1 SN (100 yr)~ 1 (10 10 L^T 1 ) given in 
Table 4 in Cappellaro et al. (1999). The B-band luminosity of the 
Milky Way is ~ a few times 10 I0 L Q and therefore our calibration 
seems to be about right. 

5.3 Comparison with other studies 

The expected LIGO/VIRGO detection rates can be determined ei- 
ther from binary population synthesis calculations, or from obser- 
vations of Galactic NSNS systems (binary pulsars). Both methods 
involve a large number of uncertainties - see the excellent review 
by Kalogera et al. (2001) and references therein. Current estimates 
for the Galactic merger rate of NSNS systems converge between 
10~ 6 - 10~ 4 yr _I . Our results belong to the lowest end of this inter- 
val. However, notice from Table that if (when) other population 
synthesis codes start using real /1-values (rather than a constant) 
their merger rates will go down by a factor of ~10. 

In general there is no consensus on the formation ratios of 
the various compact binaries. This is due to the /1-values discussed 
in this paper, the helium star evolution and the unsettled question 
of hypercritical accretion. Studies including hypercritical accretion 
adopt different approaches and do not agree on the extend and out- 
come of this effect. Studies including this process turn many of the 
(first-born) neutron stars into black holes, resulting in a much larger 
ratio of BHNS/NSNS. In our study we find a relatively high num- 
ber of BHBH systems being formed. We wonder if other codes re- 
ject systems which do not evolve through a (semi-detached) contact 
phase (i.e. if the primary star evolves to giant dimensions without 
filling its Roche-lobe). As we have seen in this study, the direct- 
supernova mechanism (Kalogera 1998) forms a significant fraction 
(~l/3) of all close-orbit BHNS systems and also contributes (~5%) 
somewhat to the formation of NSNS mergers. 

The distribution of separations (or lifetimes) of the merg- 
ing binaries is roughly in agreement with the results of other re- 
cent studies. The double neutron stars in the study of Belczynski, 
Kalogera & Bulik (2002) have slightly tighter orbits. This is mainly 
due to their higher limit of M^" 1 resulting in more systems going 
through a second common envelope phase. Their BHNS systems 
are marginally wider than ours, since their progenitors all are more 
massive than this limit. The distribution of NSNS separations found 
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by Bloom, Sigurdsson & Pols (1999) is close to ours, even though 
they have not included mass transfer from a helium star. This is 
somewhat a surprise since it is required that the helium stars lose 
their envelope in order to form very tight orbits without merging. 

We find our distribution of space velocities of the merging 
binaries to closely resemble both that of Belczynski, Kalogera & 
Bulik (2002), and that of Bloom, Sigurdsson & Pols (1999). Our 
distribution of merger sites mostly rely on the lifetimes and space 
velocities, and is therefore also in agreement with their studies. 

Our parameter variation studies show that the distributions of 
lifetimes and space velocities of the merging binaries are quite sta- 
ble (unlike the formation rates). This might seem somewhat sur- 
prising. If, as an example, the inspiral of the neutron star during 
the common envelope is more efficient, then it would seem that the 
distribution of NSNS systems should become tighter. This does not 
happen due to the fact that only a limited range of post-inspiral dis- 
tances lead to survival of NSNS systems. If the orbit becomes too 
tight the neutron star merges with the core of the giant. If, on the 
other hand, the orbit is too wide then the supernova explosion is 
likely to disrupt the binary. 



5.4 Progenitors of short-duration GRBs 

As discussed in the introduction, it seems likely that one needs a 
disk of nuclear matter around a black hole in order to produce a 
GRB. Such a disk can be formed by the collapse of the core of a 
massive star or by a merger of either a double neutron star or a 
mixed neutron star/black hole system. The merger of a black hole 
and a white dwarf (Fryer et al. 1999), or the merger of a neutron star 
or black hole with a CO- or helium core (Fryer & Woosley 1998), 
as a result of spiral-in in a CE, may not be a good candidate for pro- 
ducing a short-duration GRB. One concern is that a white dwarf, or 
a CO/He-core, is a relatively low-density object compared to nu- 
clear matter and it has a relatively large dimension (> 5000 km). 
The resultant disk would probably then be accreted much slower 
than needed for making a short-duration GRB. Merging double 
black holes are not candidates either due to the complete lack of 
any accretion disk. Therefore, we only considered the merging of 
double neutron stars (NSNS) or neutron star/black holes (BHNS 
and NSBH) as progenitors of the short-duration GRBs. 



5.5 Do all GRBs origin from binaries ? 

It is interesting to notice that most fireball models eject only a little 
matter (~ 10~ 5 M Q ) and require a rotating progenitor. To minimize 
the amount of ejected matter, collapsars require that the hydrogen 
envelope of the massive star has been removed before the collapse 
(MacFadyen & Woosley 1999). So the progenitors are in fact ro- 
tating helium stars (Wolf-Rayet stars) collapsing into black holes. 
Strong winds could remove the hydrogen envelope but may also 
remove too much spin angular momentum. An effective way of re- 
moving the envelope is via RLO in a binary system. Hence, it is 
possible that all GRBs origin from massive interacting binaries. In 
a recent paper (Lee, Brown & Wijers 2002) it has been suggested 
that indeed collapsars (hypernovae) are formed in tight binaries and 
that soft X-ray transients (SXTs) with black hole accretors are the 
descendants of GRBs. 



5.6 Hypercritical accretion 

It has been argued that a neutron star engulfed in a common enve- 
lope might experience hypercritical accretion and thereby collapse 
into a black hole (e.g. Chevalier 1993; Brown 1995). However, this 
idea seems difficult to conciliate with the observations of a num- 
ber of very tight-orbit binary pulsars (Tauris, van den Heuvel & 
Savonije 2000). For example, the system PSR J1756-5322 which 
has a CO white dwarf companion (~ 0.7 M Q ) with an orbital pe- 
riod of only 0.45 days. From an evolutionary point of view there 
is currently no other way to produce such a system besides from 
a CE and spiral-in phase. From a theoretical point of view it has 
been argued that the hypercritical accretion can be inhibited by ro- 
tation (Chevalier 1996) and strong outflows from the accretion disk 
(Armitage & Livio 2000). A number of other population synthe- 
sis studies have included hypercritical accretion as always leading 
to collapse of the neutron star (e.g. Fryer, Woosley & Hartmann 
1999), or as an intermediate approach in which collapse only hap- 
pens if the neutron star accretes a certain amount of matter (e.g. 
Belczynski, Bulik & Kalogera 2002). In the latter case a new pa- 
rameter, the critical mass for a neutron star to collapse into a black 
hole is introduced (depending on the equation of state of neutron 
stars). Note however, that the masses of neutron stars determined 
in all of the five detected double neutron star systems are quite sim- 
ilar and close to a value of 1.4M Q . Hence, in all these cases it is 
clear that the first-born neutron star did not accrete any significant 
amount of material from the common envelope. For the reasons 
mentioned above, we currently refrain from the possibiliby of hy- 
percritical accretion onto neutron stars which successfully eject the 
envelope of their companion star. 



6 CONCLUSIONS 

We have performed binary population synthesis in order to esti- 
mate the properties of massive compact binaries. In particular we 
have focused on the kinematics of close-orbit binaries which merge 
within a Hubble-time as a result gravitational wave radiation. Our 
simulated populations of massive double degenerate systems are 
relatively stable against variation of the parameters of the code. 

On a number of points we find that our description is more 
advanced than the work previously published in the literature: 

• We have used realistic values of the /1-parameter of the 
common envelope and spiral-in evolution. We estimated this 
parameter for each individual binary donor star. 

• We have used realistic models for the evolution of helium stars 
and included the effects of mass transfer from helium stars. 

• We have included the effects of energy feedback from an 
accreting compact object in the common envelope evolution. 

• We have applied reduced kicks to black holes in accordance 
with observations of black hole binaries. 

We find that the formation rate and merger rate of massive 
double degenerate binaries is dominated by double black holes sys- 
tems. The merger rates and the number of such systems present in 
our Galaxy today are estimated to be (based on model A): 

• The Galactic merger rate of BHBH systems is ~ 9.7xl0~ 6 yr~' 
and the present disk population is ~ 320 000 systems 

• The Galactic merger rate of BHNS systems is ~ 5.0xl0~ 7 yr~' 
and the present disk population is ~ 6 900 systems 

• The Galactic merger rate of NSBH systems is ~ 8.4x 10~ 8 yr 1 
and the present disk population is ~ 6 400 systems 
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• The Galactic merger rate of NSNS systems is ~ 1.5x10 6 yr 1 
and the present disk population is ~ 3 700 systems 

The reason for the relatively low value of e.g. the NSNS merger 
rate is a result of our use of realistic /l-values. The mixed neutron 
star/black hole binaries have somewhat lower formation rates than 
the NSNS binaries, but a slightly larger population present in the 
Milky Way today, as a result of their longer merging timescales. 
Our studies also reveal an accumulating numerous population of 
very wide-orbit BHBH systems which never merge. 

We estimate the detection rate for LIGO II and find that this 
rate is totally dominated (> 99 %) by BHBH mergers: 

• -840 yr~' for the BHBH mergers 

• ~3.4 yr~' for the BHNS mergers 

• ~0.6 yr~' for the NSBH mergers 

• -2.0 yr~' for the NSNS mergers 

With a bit of luck LIGO I may detect a merging event of a BHBH 
system (0.25 yr -1 ). The numbers quoted above are from model A. 
For other reasonable choices of parameter values these numbers 
typically differ within a factor of ~ 2 (see Table0. 

We have also explored the movement of the binaries in the 
gravitational potential of galaxies. This allows us to estimate the 
distribution of merging offsets from the centers of the galaxies. We 
find that a large fraction of the compact binaries are retained within 
the massive galaxies, but for lower galaxy masses these binaries 
are able to escape and merge outside the galaxies. This result is 
not consistent with the observed offset distribution of long-duration 
GRB afterglows. We therefore conclude that merging compact ob- 
jects cannot account for the ensemble of GRB offset distributions 
observed today - e.g. simulated NSNS binaries have a median pro- 
jected radial (offset) distance of ~ 4 kpc with respect to their host 
galaxy, while the observed value from 20 long-duration GRBs is 
only ~ 1.3 kpc. 
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